arXiv:1507.08860v2 [physics.bio-ph] 17 Mar 2016 


The influence of dispersal on a predator-prey system 

with two habitats 
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Abstract 

Dispersal between different habitats inflnences the dynamics and stability of 
populations considerably. Furthermore, these effects depend on the local in¬ 
teractions of a population with other species. Here, we perform a general and 
comprehensive study of the simplest possible system that includes dispersal 
and local interactions, namely a 2-patch 2-species system. We evaluate the 
impact of dispersal on stability and on the occurrence of bifurcations, includ¬ 
ing pattern forming bifurcations that lead to spatial heterogeneity, in 19 dif¬ 
ferent classes of models with the help of the generalized modelling approach. 
We hnd that dispersal often destabilizes equilibria, but it can stabilize them 
if it increases population losses. If dispersal is nonrandom, i.e. if emigration 
or immigration rates depend on population densities, the correlation of sta¬ 
bility with dispersal rates is positive in part of the models. We also find that 
many systems show all four types of bifurcations and that antisynchronous 
oscillations occur mostly with nonrandom dispersal. 
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1. Introduction 


In ecology both the exploration of dynamical models of food webs (Pas- 


cual and Dunne| 

2005 

Thompson et ah, 

2012; 

Rooney and McGann 

2012 
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and the study of spatial metapopulations(Holland and Hastings, 2008; Han 


ski and Gaggiotti, 

2004 

Tromeur et ah 

2013) 


2013) are well established lines of 


research. The two modelling approaches emphasize different aspects of eco¬ 
logical dynamics that are both relevant for most species: the dispersal be¬ 
tween different habitat patches and tropic interactions with other species. 
Yet, models that combine dispersal with trophic dynamics have only recently 
begun to appear. In the following we refer to such models that include both 
these features as meta-foodwebs. 

An elegant meta-foodweb model that is based on rates for colonization 
and extinction was proposed by Pillai (Pillai et ah, 2011). More detailed 


dynamical models were proposed in Abrams and Ruokolainen (2011), Abdl- 
laoui et ah (2007) and Jansen (2001|), to name a few. The different models 


draw motivation from different biological systems and hence make different 
modelling choices regarding the nature of dispersal and local dynamics. For 
example Jansen (2001) assumes diffusive dispersal between patches which 
is appropriate for simple life forms such as bacteria, whereas [Abrams and 


Ruokolainen (2011) make dispersal dependent on growth rate differences be¬ 


tween patches, implying that individuals can actively choose the site with 
the best growth conditions. Though many of these studies only describe 
two patches, conclusions from multi-patch models are often consistent with 
those from two-patch models (Jansen and de Roos, 2000), and therefore the 
insights gained from 2-patch systems have wider applications. 

The two mentioned examples for implementing dispersal are only a small 
subset of the large space of possibilities. In the literature different assump¬ 
tions are made regarding the functional forms of the number of emigrants 
from a given habitat patch, the choice of destination, the proportion of sur¬ 
vivors that arrive as immigrants in the destination patch and the settle¬ 


ment success (Amarasekare, 2008 Armsworth and Roughgarden, 2008; Row¬ 


ell, 2010). In the simplest case a fixed proportion of the population emigrates 


per unit time and instantaneously and losslessly settles in a randomly chosen 

This type of migration is usu- 
and often leads to 


2004). 


neighbouring patch (Leibold et ah 
ally called “random dispersal” or “diffusive migration” 
a synchronisation of the population dynamics of the two patches (Goldwyn 
and Hastings, 2008; Jansen| 2001). Examples for non-random dispersal are 


predator evasion and predator pursuit (Li et ah, 2005), or a migration rate 


that is proportional to the difference in growth rates between two patches 
(Abrams and Ruokolainen, 2011). 

The type of dispersal strongly affects the stability and the dynamics of the 
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system. In general, more rapid dispersal is more likely to synchronize popula¬ 
tions, although synchronisation does not necessarily require strong dispersal 
(Liebhold et ah, 2004). With adaptive dispersal, antisynchronous oscillations 
of the two patches are observed, which increase metapopulation persistence. 
The less similar populations are the less likely is it that dispersal is synchro¬ 
nizing (Esa Ranta and Lundberg, 1998). Other authors hnd that increased 


dispersal can decrease synchrony in population dynamics depending on the 
interactions between migrating species (Koelle and Vandermeer, 2005). An 
general investigation of metapopulations based on a linear stability analysis 
(Tromeur et ah, 2013) found that costly dispersal and social fencing are sta¬ 
bilizing, while positive density dependence and settlement facilitation reduce 
stability. Other papers have shown that costly dispersal might be destabi¬ 


lizing to a metapopulation with homogeneous patches (Kisdi, 2010) so the 
specihc mechanisms appear to be of importance. The effects of dispersal on 
stability can depend not only on the type but also on the intensity of the 


dispersal (Briggs and Hoopes, 2004) 


While much progress has been made for metapopulations and for specihc 
example systems for meta-foodwebs, a broad and general understanding of 
how different factors impact the stability of meta-foodwebs is still lacking. 
For instance is it unclear under which conditions dispersal has a stabilizing 
impact. Furthermore, meta-foodwebs can potentially undergo various types 
of instabilities. The study of foodwebs has provided abundant examples of 
two basic mechanisms of instability. The saddle-node bifurcation, which can 
lead to the relatively sudden collapse of populations, and the Hopf bifurca¬ 
tion, which gives rise to (at least transient) oscillations. In meta-foodwebs 
both of these instabilities can occur in two variants. The hrst of these af¬ 
fects all patches equally and is thus closely related to the bifurcations in 
non-spatial food web models. In the second type different patches are af¬ 
fected differently. They are thus reminiscent of pattern-forming instabilities 
such as the Turing and wave-instabilities, which are known from systems of 
partial differential equations ( Segel and Jacks^ 1972). While also these 
bifurcations lead to instability, their impact on the overall population den¬ 
sity is less pronounced, and they act as drivers of heterogeneity, which, in 
the long run, might beneht the system. In addition to these four types of 
instabilities, meta-foodwebs show further instabilities that occur out of the 
attractors created by these basic bifurcations, such as bifurcations involving 
heterogeneous hxed points, and nonlocal bifurcations involving limit cycles 
or strange attractors. 
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In order to gain an overview of the possible dynamical patterns of a sys¬ 
tem and their requirements, the generalized modelling approach 

), which is based 
powerful. The idea behind 
this approach is to consider models where the kinetics of some processes 
have not been restricted to specihc functional forms. Considering a model 
with a specihc structure, but containing general functions, allows capturing 
well-known insights into the structure of the system, without requiring often 
questionable assumptions on the exact form of kinetics. Further advantages 
are a short computation time and ease of biological interpretation. We will 
conhne our study to bifurcations out of homogeneous equilibria. This means 
that instabilities of heterogeneous systems and nonlocal bifurcations are not 
considered. 

In this paper, we investigate the dynamics of two species on two identical 
patches using the generalized modelling approach. We are able to analyse 
a broad class of models that includes several previously studied systems as 
special cases. We focus on the effect of the type and strength of migration on 
the stability and the dynamics of the system. We hnd migration in most cases 
to be either destabilizing or to have a marginal effect on stability. However, 
complex migration rules allow for a stabilizing inhuence of dispersal and 
can produce saddle-node and Hopf bifurcations and spatial-pattern forming 
bifurcations. 

2. Model 

2.1. Generalized modelling formulation 

We consider a system consisting of two habitat patches, where each patch 
i can potentially sustain a prey population Xi and a predator population Fj. 
We assume a homogeneous system, such that both patches are described by 


Feudel, 2006t Yeakel et ah, 2011 Gross et ah, 2009 


stability analysis of steady states, is particularly 


(Gross and 
on a linear 
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identical parameter values. The population dynamics are described by 


Xi = G{X,)-K{X,)-F{XuY,) 

+7]^E^{X2, F2, Xi, Fi) - Fi, X 2 , F 2 ) 

F2 = AF(Xi,Fi)-D(Fi) 

+r]^E^iX 2 , F2, Xi, Fi) - E^iX,,Y,, X2, F2) ( 1 ) 

X 2 = G'(X2)-X(X2)-F(X2,F2) 

+r]^E^{Xi, Fi, X2, F2) - E^{X 2 , F2, Xi, Fi) 

F 2 = AF(X2,F2)-D(F2) 

+r]^E^{X,, Fi, X2, F2) - E^{X 2 , F2, Xi, Fi), 


where we used the dot over a variable to indicate the temporal derivative. 
The variables are in arbitrary units. For the purpose of this paper we will 
assume that they describe the system in terms of carbon biomass density, 
however, the same equations also apply to other measures of population, such 
as abundance. The prey population density changes due to a growth rate 
G{Xi), a respiration/mortality rate K{Xi), and a rate of biomass loss by 
predation E{Xi,Yi). Predator populations have a growth term XE{Xi,Yi), 
with the prefactor A describing the efficiency of the energy conversion. The 
respiration/mortality rate of the predator is given by D{Yi). The rate of em¬ 
igration is E^(K,Y) for both species U = X,Y and the migration loss factor 
is r]^. In the most general case, emigration rates depend on all four popula¬ 
tions. The case where the emigration rate of population Ui is proportional 
to Ei and independent of other variables corresponds to diffusive migration, 
otherwise we get different versions of adaptive migration. 

Modes of the form Eqs. Q can have multiple feasible steady states, de¬ 
pending one the choice of functional forms and parameter values. In the 
generalized model we cannot compute the steady states. However, a central 
insight is that we can still compute conditions for the stability of steady 
states and express them in the form of meaningful ecological parameters. 
For this purpose we consider an arbitrary feasible, but not necessarily stable 
steady state. 

The normalized biomasses are 

Xi Fi 

~ Y*' 


( 2 ) 


5 


and the normalized fnnctions are 


H(X*,Y*) H* 


(3) 


with H = E,F,G,K and the asterisk (*) is nsed to denote the valnes of 
variables and fnnctions in this steady state nnder consideration. 

In terms of these normalized qnantities, Eqs. ([^ take the form (with 
fi = {xi,yi) and f 2 = ( 0 : 2 , 1 / 2 )) 


G 


K* 


EX*^x 

X* 


Xi 

2/1 

X2 


- —f{xi,yi) + 


e^(r2,fi) - ^^e^(fi,r2) 


E 


X* 


X, 


X* 


?Y*^Y 




e^{r2,fi) 


E 




Y* 


Y* zy* Z7* 

^9ix2) - ^Hx 2) - -^f{x2,y2) + 


)^p* JJ* TPY*^Y 


ypfix2,y2) - y^d{y2) + 


Y=^ 


E^*V^ „x 
X* 

E 


e {ri,f 2 ) 
E 


(4) 




e (ri,f 2 )-^—e (f 2 ,ri) 


X* 


y* 


e (nX 2 ) - y;ye (r 2 ,fi) 


It is now nsefnl to identify the total biomass tnrnover rate of the popnlations 
at their steady state. At the steady state, all the fnnctions in (|^ take the 
valne 1, and the gain and loss terms are eqnal and can be denoted as 


a 


X 


a 


Y 


G* E^*y^ 

-^- — 

X* X* 

XE* E^*y^ 
-^- — 


D* E^* 

Y* Y* ' 


(5) 

( 6 ) 


If the variables describe abnndances then the parameter is the total 
tnrnover rate for the prey in a patch. For instance, a valne of 0.25/year 
wonld indicate that an individnal spends on average 4 years in the patch 
before either dying or emigrating. 

It is convenient to measnre time in terms of the inverse of the prey 
tnrnover rate . In these rescaled time nnits the tnrnover rate of the prey 
is 1 and the tnrnover rate of the predator is 


a 


Y 


a = 


a 


X 


(7) 


Fnrthermore, we denote the different relative contribntions to the growth and 
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loss terms by so-called scale parameters 5, u and p and their complements 
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v 


u 


p 


u 


1 F* 



p^ X* ' 

1 

EU*^u 


U* 

1 

EU* 


u* ’ 


6={l-S) 


1 K* 


= (1 - 


1 {G*-,XF*} 
^ U* 


p^ = (l-p^) = 


U\ 


1 {K* + F*;D*} 

U* 


( 8 ) 

(9) 

( 10 ) 


Scale parameters describe the branching of the biomass flow, i.e., they 
are proportions of the total biomass inflnx or ontpnt attribnted to a specihc 
fnnction or mechanism. The parameter 6 denotes the proportion of energy 
intake of the prey that is eventnally lost again dne to respiration or mortality, 
whereas the 6 is the proportion of the energy intake that is eventnally lost dne 
to predation. The parameters u and 9 denote the relative contribntions of 
migration and feeding to the total gain, respectively. The parameters p and 
p are the connterpart to v and 9 and denote the relative loss by emigration 
and by the within-patch processes (respiration/mortality and predation). A 
valne of zero of a scale parameter means no biomass flow attribnted to a 
mechanism (e.g. 5=0 means no gain by predation) and valne of 1 means a 
biomass flow completely dominated by a mechanism (e.g. p^=l means loss 
only by predation and respiration). 

Using the tnrnover rates and the scale parameters, we can write Eq. (|^ 
as 


xi = a^[9^g{xi) - pfSk{xi) - p^6f{xi,yi) + u^e^{^,y) - p^e^{^,y)] 
yi = a^[9^ f{xi, yi) - pjd{yi) + z/^e^(x,y) - p^e^(x,y)] (11) 

±2 = a^[9^g{x2) - p^5k{x2) - p^5f\x2,y2)+ - p^e^{-K,y)] 

y 2 = a^[9^f{x2, 1 / 2 ) - P^d{y 2 ) + z/^e^(x,y) - p^e^(x,y)]. 


For analysing the stability we linearise the system aronnd the steady state 
nnder consideration. In the normalized system the steady state is at x=y=l. 
The linearisation can then be expressed in terms of the Jacobian matrix. For 
a system of fonr dynamical variables this is a 4 x 4 matrix dehned by 


J, 




dVi 

dV, 


where V = {xi,yi,X 2 ,y 2 ) is the set of state variables, and 
the derivatives are evalnated in the steady state (1,1,1,1). 


( 12 ) 

indicates that 
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For the case of a system with two identical patches the 4x4 matrix can 
be written as a block matrix 


J = 


L M 
M L 


(13) 


with a matrix L that describes the in-patch dynamics, 


L = 


1 0 
0 a 


— p^u' 

~Y I Y ''Y Y Y 
U 'y + U K — p K 


-p^6lp H- — p^ 

(14) 


~Y I ~Y Y I Y ''Y Y Y 
V Ip — p p + V OJ — p UJ 


and a matrix M that captures between-patch dynamics, 
M = 


1 0 
0 a 


j,X^X_pX^X ^X^X_pXf^X 

U^K^-P^k^ U^U^-P^U^ 


(15) 


The new parameters in these matrices are called exponent parameters, and 
they are the derivatives of the population dynamic functions with respect to 
the prey or predator populations. The dehnitions of the exponent parameters 
are 


dg{x) 


7 : = 


dx 
df{x,y) 


P := 


dx 


tp : = 


dx 
df{x,y) 


Y __ dd{y) 

jl 


dy 


u^: = 


k^: = 


u^: = 


k^: = 


de^ixi,yi,X2,y2) 


dx2 

de^ixi,yi,X2,y2) 


dy2 

de^{xi,yi,X2,y2) 


u^: = 


:= 


dy 


de^{xi,yi,X2,y2) 


dxi 

de^ixi,yi,X2,y2) 


(16) 


dy2 

de^{xi,yi,X2,y2) 


dxo 


u^: = 


K^: = 


dyi 

de^{xi,yi,X2,y2) 


dyi 

de^{xi,yi,X2,y2) 


dxi 


The exponent parameters are so-called elasticities, i.e. they are logarithmic 
derivatives of the original functions. For example 


dg{x) 


dx 


d\og{G{X)) 


01og(Jf) 


( 17 ) 



Such logarithmic derivatives have a number of advantageous properties. Elas¬ 
ticities as a measure of nonlinearity were hrst introduced in the 1920s in eco¬ 
nomic theory, because of their statistical properties which allows them to be 
estimated precisely based on limited noisy data. For the same reason these 
parameters are now commonly used in metabolic control theory. 

In generalized models we use elasticities mainly because they allow for 
an intuitive interpretation: For any linear function, say G'(X) = AX, the 
corresponding elasticity (91og(G(X))/(91og(X)|i is 1, regardless of the slope 
A. More generally, for any power law G/{X) = AX^ the elasticity is the 
power law exponent p. For more complex fnnctions the elasticity provides an 
intuitive measure of the degree of saturation of the nnderlying process. For 
example for the Rolling type-II kinetics the corresponding elasticity is close 
to 1 in the linear regime at low prey density, bnt approaches zero at high 
prey density, where predators are satnrated. 


Parameter 

Interpretation 

Range 

Examples 

0 

Elasticity of the gain 
fnnction 

[0,1] 

0: growth independent of pop¬ 
nlation size, e.g. dne to nntri- 
ent limitation, 1: growth prop, 
to population density 

/i^ 

Sensitivity of prey 
mortality to prey 
popnlation 

[1,2] 

1: constant mortality 2: mor¬ 
tality proportional to density, 
e.g. dne to diseases 


Sensitivity of predator 
mortality to predator 
population 

[1,2] 

same as for prey 

7 

Sensitivity of preda¬ 
tion gain to prey pop¬ 
ulation 

[0,2] 

Rolling type fnnctions: value 
close to 1 for small prey popu¬ 
lation, low valne for large pop. 
dne to satnration 


Sensitivity of preda¬ 
tion gain to predator 
population 

[0,1] 

1: No predator competition, 
lower valnes are dne to preda¬ 
tor interference 


Table 1: The five exponent parameters that characterize in-patch dynamics, their meaning, 
and their typical range. 


The uj and k are the migration exponent parameters, where the exponent 
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bj describes the dependence of the emigration rate on the density of the 
emigrating population, whereas k captures the dependence of emigration on 
the density of the other species in the patch, respectively. The parameters 
with a hat a), k are dehned analogously but describe the dependence on the 
densities in the destination patch. Different types of adaptive migration 
imply different parameter ranges for these exponent parameters. Simple 
diffusive migration for species U means = 1, and when the emigration 
rate increases with population size we obtain an exponent uP > 1. In the case 
of predator evasion, prey emigration increases with predator density, which 
means a value > 0. In predator pursuit, predator emigration decreases 
with prey density, which means < 0. 

The remaining hve parameters (hrst hve in Eq. (16)) describe the elas¬ 
ticity of the local (within-patch) processes. These parameters have been dis¬ 
cussed extensively in previous work (e.g. Gross and Feudel (2006)). Hence, 
we summarize their interpretation in Tab. 


2.2. Linear stability, eigenvalues, and bifurcations 

In order to evaluate the stability of a steady state, we have to calculate 
the eigenvalues of the Jacobian at the steady state. Due to the symmetric 
block structure of the Jacobian, the eigenvalue equation can be solved with 
the ansatz(MacArthur et al., 2008) 


^ Li L2 Ml M2 ^ 


( ^ 


( ^ 

L3 L4 M3 M3 


6 

= A 


M\ M2 L\ L2 


±6 


±6 

^ M3 M4 L3 L4 y 


\ ±^2 / 


\ ±^2 / 


with Li and Mi denoting the matrix elements of L and M. This is equivalent 
to solving the 2x2 problem 

(L^±M^ ± m 2 \ ^ / eM 

\Li,±M^ L,±M, ) ' 


Every solution of the eigenvalue equation describes an eigenmode of the sys¬ 
tem, i.e., a specihc perturbation that retains its shape while it grows or 
declines in time. For the plus sign, the eigenvector components relating to 
the two patches are identical, for the minus sign, they point in opposite 
directions. 
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The steady state is stable if the real parts of all eigenvalues are negative. 
The eigenvalues for a 2x2 matrix J' can be written in terms of trace (T) and 
the determinant (A) of the matrix as 


A 



T(J') ± 



A(J'). 


( 18 ) 


In order to obtain a stable steady state, the real part of both eigenvalues must 
be negative, which is the case when the trace is negative and the determinant 
is positive. These criteria must be satished for both matrices L ± M (below 
denoted by the index ”+” for ”L+M” and for ”L-M”). 

As parameters are changed, a stable steady state can become unstable when 
the parameter change causes an eigenvalue to cross the imaginary axis and 
acquire a positive real part. For our system, there are four different types of 
such local bifurcations. 

Let us hrst consider the ” J+ = L+M” case. Instabilities that are detected 
by the analysis of this matrix affect both patches equally and in synchrony. 
First, stability of the steady state can be lost in a saddle-node bifurcation, 
which occurs when A(J_|_) = 0 and T(J_|_) < 0 while both eigenvalues of J_ 
have a negative real part. In this bifurcation a sudden change in the pop¬ 
ulation densities occurs, which most likely results in the collapse of one or 
both populations in both patches simultaneously. The second fundamental 
bifurcation in which stability can be lost is the Hopf bifurcation. This bi¬ 
furcation occurs when T(J_|_) = 0 and A(J+) > 0 while the eigenvalues of 
J_ have a negative real part. The Hopf bifurcations detected in the L-l-M 
matrix gives rise to synchronous oscillations. From the study on non-spatial 
predator-prey systems it is well known that oscillation amplitudes can grow 


rapidly after the bifurcation, leading to subsequent extinctions (Rosenzweig 


and MacArthur, 1963). 


Let us now consider the J_ = L — M matrix. Bifurcations detected in 
this matrix affect the two patches in opposite directions. Again, stability 
can be lost either in a saddle-node or in a Hopf bifurcation. However, now 
the Saddle-node bifurcation leads to a shift where each population increases 
in one patch and decreases in the other. This occurs when A(J_) = 0 and 
T(J_) < 0 and both eigenvalues of J+ have a negative real part. A Hopf 
bifurcation detected in the L — M matrix leads to anti-synchronous oscilla¬ 
tions. While this can lead to large-amplitude oscillations in both individual 
patches, the overall biomass in the system will stay nearly constant as the 
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loss in one patch is compensated by gains in the other. This type of bifurca¬ 
tion occurs at T(J_) = 0 and A(J_) > 0 with both eigenvalues of J+ having 
a negative real part. 

Both bifurcations occurring in the L — M matrix can be considered as 
pattern-forming bifurcations that create spatial heterogeneity. In fact, the 
saddle-node bifurcation in the L — M matrix is closely reminiscent of the Tur¬ 
ing bifurcation in systems of partial differential equations, which gives rise 
to stationary patterns. The Hopf bifurcation in the L — M matrix is rem¬ 
iniscent of the wave instability bifurcation in partial differential equations, 
which leads to travelling waves. 

In the following we will use the terms Hopf and saddle-node bifurcation 
only for the corresponding bifurcations from the L+M matrix. For simplicity 
we will denote the respective bifurcations in the L — M system as Turing and 
wave instability. We emphasize that this is a loose usage of the bifurcation 
names that we adopt here as it leads to the right ecological intuition although 
it is not strictly mathematically justihecQ 

2.3. Classes of models studied in the following 

In order to evaluate the proportion of stable systems and the frequency 
of the different types of bifurcations, we need to specify intervals for the 
exponent parameters and the scale parameters. Different classes of models 
are characterized by different choices for these intervals. We use several 
models from the existing literature as well as more general models. All these 
models are listed in Tab. [U 


^Strictly, a bifurcation is of a given type only if it can be reduced to the type’s normal 
form. For instance the saddle-node bifurcation observed in the L — M system can be 
reduced to the saddle-node normal form but not to the Turing normal form. So strictly 
it is a saddle-node and not a Turing bifurcation. Still one can justify denoting this bifur¬ 
cation as ’Turing’ as follows: The bifurcation would not change fundamentally when we 
considered a system with more than two patches. The bifurcation would then be governed 
by a matrix that is closely reminiscent of the network laplacian, which in turn can be seen 
as a discretization of the real space laplacian on a complex network. In this sense even the 
2-patch system can be interpreted as a discretization of an underlying continuous space in 
which the bifurcation would be a true Turing bifurcation. 
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Group (1) represents the most general models, for which all exponent 
parameters have their maximum meaningful range. In group (2), migration 
is diffusive for both species, but the other local exponent parameters are still 
the in maximum meaningful range. Models (4) to (7) have a Holling type 
2 functional response with logistic growth and no predator interference, and 
they differ with respect to the migration term. Model (4) Jans95 shares the 


local patch dynamics with most of the models in the literature (Rosenzweig 


and MacArthur| 1963) and has a simple diffusive migration. Model (5) is the 


same as model (4) but with no prey migration. Model (6R) Huang adds that 
over-abundance of prey facilitates predator migration while model (6) incites 
predators to stay if there is plenty of food available. Model (7) Abramsll 
implements predator pursuit, which means that predators migrate towards 
the patch with larger prey abundance. 

Models (9) and (10) establish a relation between consumption rate and 
migration: the migration rate of the predator scales in the same way as its 
feeding rate (i.e., the biomass production), and the migration rate of the prey 
scales with its own growth rate as well as with the predator feeding rate and 
with predation on the other patch. 

We group all these models into two classes, with class I comprising the 
models with general intervals for the hve in-patch scale parameters (models 
(1), (2), (9)) and class II comprising models (4) to (7) and (10), which are 
based on the Rosenzweig-MacArthur model. 

Models (3) Mchich and (8) ElAbdllaoui £x more parameters for the in¬ 
patch dynamics but have interesting migration rules: Mchich lets the prey 
flee if there are many predators in the own patch. ElAbdllaoui additionally 
makes the predator sedentary if there is plenty of food available. 


3. Results 

3.1. Proportion of stable systems 

We evaluated the proportion of stable states in an ensemble where the 
scale and exponent parameters were chosen at random from the intervals 
indicated in Tab. We generated 10^ random sets of parameters for each 
model where each parameter was drawn uniformly from the respective in¬ 
tervals. We define the proportion of stable webs (PSW) as the number of 
parameter sets which produce a stable steady state, divided by the total num¬ 
ber of sets sampled. We use a random distribution of values for exponent 
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and scale parameters as we do not assnme any relationship between parame¬ 
ters. A specific scenario will have most likely dependencies between fnnctions 
(e.g. predator feeding rate and prey mortality share common variables) bnt 
these are snbsets of onr parameter space. Weights for the parameter dis- 
tribntion wonld skew the results towards particular assumptions which we 
want to avoid. In the absence of migration (i.e., for single patches), we ob¬ 
tained PSW=0.939 for class I, PSW=0.639 for class II and PSW=0 for model 
Mchich and PSW=1 for model ElAbdllaoui. These results are a measure of 
the stability of the within-patch dynamics alone. We know from Levins 


(1974); Gross et ah (2009) that mortality and feeding terms tend to be pos¬ 
itively correlated with stability and growth terms are negatively correlated 
with stability for local dynamics. This means that high values for 7 and /i*’^ 
and low values for cj) and ip enhance stability. Based on this knowledge we 
expect model Mchich to be less stable than model ElAdbllaoui, and systems 
with class I local dynamics to be more stable than systems with class II local 
dynamics. This is confirmed by our results. Model Mchich is a special case 
as the particular parameter choices of Mchich et al. (2007) create a pair of 
purely imaginary eigenvalues, both for the case with and without migration. 
This does not allow for a conclusive linear stability analysis, and higher order 
terms are needed to judge the stability of a steady state. 

The PSW values obtained in the presence of migration are given in Tab. 

In all cases, stability with migration is smaller than or roughly equal to 
stability without migration. The proportion of stable webs ranges from 0.05 
for model (10), which has a complex migration rule, to 0.94 for class I with 
diffusive, conservative migration, i.e., models (2y) and (2z). In fact, the 
models (1), (9), (10) have the largest drop in stability due to migration. 
These models have in common that Cj and/or k are nonzero, which means that 
dispersal of a population depends on the population of the other species on 
the other patch. This can be understood by applying the qualitative stability 
considerations established by Levins (1974): Including in the Jacobian an 
element that connects different species on different patches creates a positive 
feedback loop of length 3, which has a strong destabilizing effect. 

It is interesting to note that the class II systems are less stable than the 
models with diffusive migration, ( 2 ) to ( 2 z), despite the high exponent of 
closure /ij, = 2 . A high exponent of closure is well known to be stabilizing 


since it implies a strong density limitation of the predator (Levins, 1974 Gross 


et ah, 2009 Plitzko et ah, 2012). In order to investigate the effect of the 
exponent of closure, we run the class II systems also with pa, = 1 instead of 
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2. In this case all systems became unstable. When /la- was chosen at random 
from the interval [1,2], the average stability was lower than for the isolated 
patches (around 0.50), but with trends similar to those for = 2. 


Scenario 

PSW 

PSW 


without migration 

with migration 

(1)Stand 

0.939 

0.227 

(Ix)Stand 

0.939 

0.570 

(Iz)Stand 

0.939 

0.394 

(2)Diff 

0.939 

0.924 

(2x)Diff 

0.939 

0.876 

(2y)Diff 

0.939 

0.942 

(2z)Diff 

0.939 

0.939 

(9) Growth 

0.939 

0.450 

(4)Jans95 

0.693 

0.693 

(4x)Jans95 

0.693 

0.693 

(4z)Jans95 

0.693 

0.693 

(5)JansOl 

0.693 

0.693 

(5z)JansOl 

0.693 

0.693 

(6) Huang 

0.693 

0.552 

(6R) Huang 

0.693 

0.530 

(7)Abramsll 

0.693 

0.693 

(10)Growth2 

0.693 

0.051 

(3) Mchich 

0.000 

0.00 

(3x) Mchich 

0.000 

0.00 

(8)ElAbdllaoui 

1.000 

0.625 


Table 3: The different scenarios and the proportion of stable webs with and without 
migration. The double horizontal lines separate the classes I, II, and the exceptions model 
Mchich and model ElAdbllaoui. The parameter values for the different scenarios were 
drawn uniformly from the intervals given in Tab. 

In order to obtain more detailed information about the effect of migration 
on stability, we evaluated the correlation of stability with the scale param¬ 
eters of migration. We used the Pearson product-moment correlation coef- 
hcient to compute the correlation between stability and the migration scale 
parameters. The result is shown in Fig. 
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□ Vx 

□ Vy 

°Px 

°Py 


Figure 1: The correlation of the four scale parameters of migration Px, Py) with 

stability for 10^ randomly chosen parameter sets for each model. The different models are 
arranged along the x axis, the correlation value is given along the y axis. The different 
scale parameters are coded by colour. Most models show a negative correlation of the 
migration parameters with stability. 


In most cases, the correlation is zero or negative, and there are only 
four instances of a signihcant positive correlation with migration. Positive 
correlation values mean that higher relative migration rates make the steady 
state more stable, negative correlation values indicate that the steady state 
becomes less stable. It is interesting to note that out of the hve instances 
with positive correlations four are loss terms px^y This hts together with the 
observation that loss terms are comparable to death terms, which are known 
to be stabilizing for local dynamics (Gross et al., 2009). 

To summarize, the results in Fig. [pimply that most models become less 
stable when the contribution of migration to the total gain and loss terms 
increases. However, several models become more stable with increasing mi¬ 
gration losses, ElAbdllaoui, Huang and Dijf, and only model Growth becomes 
more stable with increasing migration gains. The models ElAbdllaoui,Huang 
and Growth all show types of adaptive migration. The model Diff has simple 
diffusion as a dispersal mechanism, though both Diff and Growth have wide 
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parameter interval ranges for local parameters. 

In addition to evalnating the correlation of the scale parameters with 
stability, we also evalnated how average stability changes as a scale parameter 
is varied. Example curves that represent the qualitatively different types of 
behaviour found are shown in Fig. 

The top left graph shows a monotonous decrease, which is most often ob¬ 
tained and applies to most scenarios that have a negative correlation of the 
scale parameter with stability. The bottom right graph shows a monotonous 
increase, which is seen in the cases of positive correlation of the scale pa¬ 
rameters with stability. The bottom left graph shows an instance of a 
non-monotonous curve, found for the scenarios Stand, Growth and Growth2, 
which have complex migration rules, where more than two of the eight mi¬ 
gration exponent parameters {u , k ) are different from zero. This shows 
that cross-patch cross-species interactions as seen in these three classes can 
create positive correlation of migration scale parameters with stability within 
certain ranges even if the overall correlation remains negative. The top right 
graph is an instance where migration has no effect on stability. 
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Figure 2: Examples of the different types of functional relation between scale parameters 
and average stability. The graphs show the proportion of stable systems as function of a 
migration scale parameter. The numbers above each plot denote the model (see Tab. 
which is defined by the intervals of the exponent parameters and the relationship between 
scale parameters used to create the data. Apart from a monotonous decrease or increase, 
there are also models with a constant proportion of stable systems and those with maxima 
at intermediate values of the scale parameter. 
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Inspired by these findings, we investigated additional cases: Instead of 
linear relationships for the sensitivity of the migration on the respective pop¬ 
ulation as used in class II), we used quadratic relationships {ux,y=2) 

for the scenarios of class II. This caused clear trends, with the average sta¬ 
bility decreasing with increasing values of the scale parameters and in¬ 
creasing with px^y. Even though slope and absolute values vary, this effect 
dominates over every other influence on stability and shapes the curves al¬ 
most solely. This is understood by realizing that Ux,y = 2 implies large loss 
terms at high values of p. The dependence on v is more complicated as v 
affects only the non-diagonal elements of the Jacobian matrix. High v cre¬ 
ates larger diagonal elements for the migration sub-matrix. This increases 
the trace T(J+) and thus decreases the likelihood of a stable steady state 
even though T(J_) becomes more stable as the total stability is limited by 
the stability of each sub-matrix. 


3.2. Example of analytical computation 

In most cases, the stability curves are not easily understood. The sta¬ 
bility condition comprises four inequalities, two for each matrix LEM (see 
paragraph after eq. (18)). Each of these inequalities contains several param¬ 
eters, and checking whether they are satisfied requires the consideration of a 
multitude of cases. In the following we demonstrate for the model with the 
smallest number of free parameters, which is the model ElAbdllaouf how the 
shape of the stability curve results from the inequalities. As an example we 
use the average stability versus the scale parameter Uy, as it has a distinct 
shape of two linear sections with different slopes which can be seen in Fig. 

El 

The traces and determinants for this model are 


r(J+) = 1 - 2{1 - S){1 - px) - - px) - Px , 

T{J_) = 1-2(1-5){1-px) - 6{1-px) - px-‘2iyx-‘^ayiyy, 
A(t/+) (6( 1 -|- Px) px T ^x)(o:y(l T py ^y) ^y^y) i 

A(J_) = -2ay(l-2(1 - 5)(1 - Px) - 5(1 - Px) - Px-2.Vx)vy , 
(J( 1 -|“ Px) Px ^a;)(cij/(l T Py ^y) T • 


The steady state is stable if both traces are negative and both determinants 
are positive. The values of J, px^y and Vx,y are in the interval [0,1]. This means 
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PSW 



Figure 3: The average stability as a function of Vy for the scenario ElAbdllaoui. The solid 
red line shows the analytical result. 


that T(J+_) < 0 and A(J_) > 0. The stability is thus solely dependent on 
A(J_). Rewriting A(J+) as 


A(J+) = -ayfi ■ /s 

with fi = 6{px — 1) — px + and f 2 = l + py — 2h'y, the stability condition 
becomes /1/2 < 0, which means that /i and /2 must have opposite signs. The 
term fi is independent of Uy and thus cannot change sign as Uy is changed. 
With all parameters in fi being random numbers in [0,1], /i < 0 in 75% of 
the cases, and /i > 0 in 25% of the cases. /2 is always larger than zero as 
long as Uy < 0.5. This means that the steady state is stable in 75% of the 
cases. As Uy increases from 0.5 to 1, the probability that /2 > 0 drops linearly 
from 1 to 0. For Uy = 1, only 25% percent of the systems are stable because 
now /i must be positive. This explains the linear drop of the stability curve 
from 0.75 to 0.25. 

3.3. Bifurcations 

We evaluated how often each of the four different types of local bifurca¬ 
tions occurs as a migration scale parameter is increased from 0 to 1, averaging 
over 10^ parameter sets and over the four different migration scale param¬ 
eters. Varying the migration scale parameter means varying the relative 
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contribution of migration to the population growth and loss terms. This 
approach is different from merely increasing a migration rate in an explicit 
model and therefore leads to somewhat different statistics of bifurcations. 
The comparison with explicit models will be provided in the next section. 
The result is shown in Fig. 


Proportion of systems 


with bifurcation 



Scenario 


Figure 4: The statistics of the four types of bifurcations for all considered models for 
10® random parameter sets. The numbers on the x-axis denote the model (see Tab. [^. 
Since more than one type of bifurcation can occur as a scale parameter is increased from 
0 to 1, the maximum total height of each bar is 4. The symbols are SN for saddle-node 
bifurcation, HB for Hopf bifurcation, TI for Turing instability and WI for wave instability. 
(See Section 2.2 for the definition of the bifurcations.) 


The dominant bifurcation is always the saddle-node bifurcation. We see 
all four different types of bifurcations in models Stand and Standx as well as 
models Growth and Growth2. These are the only models in which migration 
responds to population changes differently for predators and prey. The mod¬ 
els with diffusive migration {Diff and Jans) do not show wave instabilities. 
For the scenarios of class 11 we see barely any bifurcations, with the exception 
of the models Huang,HuangR and Growth2, though model Huang only ex¬ 
hibits saddle-node and model HuangR saddle-node and Turing bifurcations. 
This is consistent with the earlier analysis of stability. If migration does not 
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cause a change in stability, there can be no bifurcations that are induced by 
migration. Comparing these results with Table we see that there is some 
correlation between the amount by which migration decreases stability and 
the number of bifurcations that occur. 

By evaluating the bifurcation statistics separately for different intervals 
of the migration scale parameters, we found that generally Hopf bifurcations 
and wave instabilities occur more often for intermediate values of the scale pa¬ 
rameters, while saddle-node bifurcations and Turing instabilities occur more 
often for large and small values. Figs. and show the corresponding his¬ 
tograms for selected models. For these hgures, we evaluated additionally the 
information whether the used parameter set would produce a stable steady 
state without migration, i.e. with all migration scale parameters set to zero. 
This subset of the total number of bifurcations is marked by the violet colour. 

Saddle-node bifurcations are more likely to occur at high migration val¬ 
ues, whereas Turing-bifurcations tend to low and mid values. The steep in¬ 
crease of the number of saddle-node bifurcations for p close to 1 is expected 
since large loss terms can lead to predator extinction, i.e., to a transcritical 
bifurcation, which is identical to a saddle-node bifurcation in the general¬ 
ized modelling approach. Synchronous Hopf bifurcations are usually found 
at mid-to-high levels of migration while wave-instabilities prefer mid-to-low 
migration. Intermediate values of the migration scale parameters imply that 
the interactions between populations within and between patches are equally 
important for the dynamics, which in turn means that the phase space in 
which the dynamical trajectories of the system evolve is four-dimensional. 
Without migration or with very fast migration the phase space is only two- 
dimensional, with less complex dynamics and less oscillations. 

It has been noted that prey pursuit or predator evasion can decrease syn¬ 
chrony (Li et ah, 2005). The scenarios ElAbdllaoui, Growth and Growth2 
implement such types of migration. Model ElAbdllaoui shows more syn¬ 
chronous Hopf bifurcations at high values of the migration scale parameters. 
Models Growth and Growth2 have no clear trend. 
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Figure 5: The number of bifurcations into oscillating states in dependence of the migration 
scale parameters, in a sample of 10® random parameter sets. The symbols above each 
graph denote the model (see Tab.[^ and type of bifurcation(”HB” stands for synchronous 
Hopf bifurcation and ”WI” stands for wave instability). The shape in the top left graph 
with most bifurcations at mid-range values is found for 38 percent of cases. Other forms 
that are seen less often (predominantly in models with complex migration dynamics) are 
shown in the other histograms. The darker bars indicate the part of parameter sets with 
a bifurcation that would have been stable without any migration. 
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Figure 6: Histograms of the number of saddle-node and Turing bifurcations over migra¬ 
tions scale parameters. The majority of bifurcations occur for scale parameters close to 
zero or 1. We sampled 10® random sets of parameters per model. The number of bifurca¬ 
tions into oscillating states in dependence of the migration scale parameters, in a sample 
of 10® random parameter sets. The symbols above each graph denote the model (see 
Tab.[^ and type of bifurcation(”SN” stands for saddle-node bifurcation and ”TI” stands 
for Turing instability). The majority of bifurcations, 67 percent of total cases, occur for 
scale parameters close to zero or 1. The darker bars indicate the part of parameter sets 
with a bifurcation that would have been stable without any migration. 
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4. Comparison with results from explicit models 

Many of the models listed in Table were investigated earlier using ex¬ 


plicit population dynamics models. These are the models by Mchich et ah 
(2007), Jansen (2001, 1995), Abdllaoui et ah| (2007), Abrams and Ruokolainen 
(2011) and Huang and Diekmann (2001). Our approach provides a different 
perspective and gives general insights. It thus complement the findings of 
models that use specific functional forms for the different gain and loss terms. 
In order to illustrate the respective strengths of the different approaches, we 
compare in the following our results with those publications. We hereby have 
to keep in mind that our study was confined to homogeneous systems with 
identical patches and to local bifurcations of steady states, while the cited 
publications often include heterogeneous systems and global bifurcations as 
well. 


Jansen (1995) sees no dependence of the bifurcation condition on migra¬ 


tion for the homogeneous steady state with both species present on both 
patches. In our study, we found (few) saddle-node and Hopf bifurcations 
in this model when the migration scale parameters are changed. Such bi¬ 


furcations are also present in the model by Jansen (1995) and are crossed 


when the parameters that characterize the in-patch dynamics are changed. 
This holds also for the follow-up paper from 2001(Jansen, 2001). This nicely 


illustrates the fact that changing one parameter in the generalized approach 
is not equivalent to changing one parameter in a compatible explicit model, 
but to changing several parameters simultaneously, and vice versa. For in¬ 
stance, increasing the migration scale parameter v in our approach means 
that a larger proportion of biomass increase is due to immigration, and that 
the biomass increase due to resource consumption decreases. In contrast, 
increasing the migration rate in an explicit model while keeping all other 
parameters fixed means that migration rate increases without a change in 
other processes. Varying scale parameters provides more generic insights as 
they refer to relative and not to absolute quantities. It is of course possible 
to change the parameters in an explicit model in a fixed relation and thus 
obtain the same kind of insight. 


Mchich et ah (2007) find a stabilizing effect of migration for a situation 


where the prey migrates with a rate that depends on predator density. This 
study uses a Lotka-Volterra model, for which a linear stability analysis is 
exact, and evaluates it in the limit of rapid migration. Our study always gives 
an eigenvalue zero of the Jacobian, since we look at homogeneous systems. 
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The Hopf bifurcations observed by Mchich et al. (2007) are due to the fact 
that they investigate the inhomogeneous case where the parameters on the 
two patches are different. The study by Abdllaoui et al. (2007) is similar 
to that by Mchich et al. (2007), but uses a type II functional response and 
hnds that migration can create limit cycles. We see a large percentage of 
stable webs, with Hopf bifurcations occurring only rarely. This leads to the 


conclusion that the majority of Hopf bifurcations seen by Abdllaoui et al. 


(2007) is due to the fact that they study an inhomogeneous system. 


Abrams and Ruokolainen (2011) state that adaptive migration produces 


frequent anti-synchronous limit cycles. Our study generally shows wave in¬ 
stabilities for scenarios with adaptive migration (models Stand, Growth), and 
such wave instabilities lead to anti-synchronous limit cycles. However, we did 
not hnd wave instabilities for model Abrams. Since lAbrams and R.uokolainenI 
(2011) do not show bifurcation diagrams, there is no contradiction between 
their and our hndings. They focussed on parameter ranges where local patch 
dynamics is oscillatory. It is well possible that in this model antisynchronous 
oscillations result only (or almost always) out of inhomogeneous hxed points 
or by non-local bifurcations. 


Huang and Diekmann (2001) investigate the case that predators cannot 


migrate while handling prey, which means that predators migrate at low prey 
densities but become immobile for large prey densities. By performing a lin¬ 
ear stability analysis, the authors found a stable steady state for a wide range 
of parameters, which can become unstable by Hopf bifurcations and saddle- 
node bifurcations. We only hnd saddle-node bifurcations when we vary the 
migration scale parameters. However, we hnd also Hopf bifurcations when we 
vary additionally one of the local parameters. This again illustrates the fact 
that changing one parameter in an explicit model corresponds to changing 
several parameters in the generalized model, and vice versa. Furthermore, 


Huang and Diekmann (2001) state that antisynchronous oscillations are al¬ 


ways unstable if they exist, in agreement with our result that a stable steady 
state cannot become unstable by a wave instability. 

To conclude this section, we would like to point out that a model with ex¬ 
plicit population dynamics explores only a small subspace of the general class 
of models that are compatible with what is known from empirical systems. 
In particular, it introduces functional dependencies between the parameters 
used in the generalized modelling approach and thus limits the space that can 
be explored. Other explicit models would lead to somewhat different func¬ 
tional relations between the coresponding generalized parameters and would 
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thus explore a different region of the space of possibilities. For instance, the 
model by Abrams and Ruokolainen (2011) does not show wave instabilities, 
which do however occur in models that belong to the same overall class of 
systems that have a growth-rate dependent migration. Also, the model by 


Mchich et ah (2007) finds Hopf bifurcations only in inhomogeneous systems. 


due to the Lotka-Volterra form of local population dynamics, while a gener¬ 
alized investigation that does not insist on this unrealistic constraint shows 
that migration can generically induce Hopf bifurcations in homogeneous sys¬ 
tems. Of course, we cannot rule out on the other hand that the generalized 
approach includes parameter combinations that cannot be satisfied by any 
realistic and explicit population dynamics model. 


5. Conclusions 

In this paper we studied the impact of dispersal on a general class of 2- 
patch 2-species predator-prey system. We used the approach of generalized 
modelling, which allows to investigate the stability of the system without 
restricting the kinetics to specific functional forms. In comparison to previous 
studies we are thus able to provide a broader overarching perspective. 

Our analysis confirms that dispersal generally decreases the local stability 
of the system. Although, dispersal creates benefits, such as the rescue effect, 
not studied here, it does not generally promote the dynamical stability of 
the system. Thus dispersal may lead to undesirable dynamics, causing for 
instance increased spatial, temporal, or spatio-temporal variability. How¬ 
ever, we found also large parameter regions in which dispersal increases the 
stability and may thus help avoiding such undesirable dynamics. 

For the case of identical patches we were able to compute the thresholds at 
which destabilization occurs analytically. Furthermore, we used a numerical 
sampling procedure for a broad survey of the impact of parameters of the 
generalized model. By restricting the parameters of the generalized model 
to appropriate ranges we were able to analyse 19 different scenarios. 

In a number of scenarios the impact of dispersal is very weak. For instance 
if dispersal occurs completely randomly it does not have a notable effect on 
stability. Moreover, density independent dispersal only affects stability if the 
growth of the prey shows effects of saturation. Finally, superlinear mortality 
rates, such as quadratic mortality, have a stabilizing effect that can be much 
stronger than the effect of dispersal, thus that dispersal does not generally 
lead to a destabilization of such systems. These findings are consistent with 








Tromeur et al. (2013), who studied the dynamics of a single population in 
a system with many patches. In contrast to this previous study we found 
that dispersal can in certain cases be stabilizing in the two-species system. 
The parameter regions in which a positive effect of dispersal on stability are 
observed are much wider when dispersal of a given species depends on the 
densities of other species in the system. 

While we studied only 2 patches, our results allow for some extrapolation 
to larger systems. We need to distinguish between pattern forming (Turing, 
wave) and non-pattern forming (Hopf, saddle-node) instabilities. In agree¬ 
ment with Abrams and Ruokolainen (2011) we find that pattern forming 
instabilities can only occur in systems with non-random dispersal, a result 
that should hold also in systems with more patches. 

|2M2 


Using general insights into the dynamics of networks (Do et al. 


MacArthur et al., 2008) we can say that non-pattern-forming instabilities 
cannot depend on network structure. If dispersal affects these instabilities 
then it does so only because it shifts the operating point of the system and in¬ 
troduces new nonlinearities, e.g. from losses during migration. Both of these 
effects can be captured faithfully in single patch models, where migration 
affects are modelled as additional terms. 

By contrast, pattern forming instabilities depend sensitively on network 
structure and thus can only be analysed if the spatial structure at hand is 
captured in the model. Previous results(Do et ah, 2012) suggest that the 


symmetric pair of patches studied here is a particularly unstable conhgura- 
tion that promotes pattern-forming instabilities, while larger, less symmetric 
systems should be tentatively more stable. 

Perhaps most importantly our hndings illustrate the complex nature of 
dispersal effects. In the simplest scenarios dispersal is neither stabilizing 
nor destabilizing, but the interplay of dispersal with other nonlinearities in 
the system can increase or reduce stability. To fully understand dispersal 
effects in complex food webs, and in particular pattern forming instabilities, 
we will eventually have to study large, many-patch, many-species systems. 
For addressing this challenge, the generalized modelling approach, used here, 
could be valuable tool. 
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